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Abstract. 

We present a finite-size scaling analysis of high-statistics Monte Carlo simulations of 
the three-dimensional randomly site-diluted and bond-diluted Ising model. The critical 
behavior of these systems is affected by slowly-decaying scaling corrections which 
make the accurate determination of their universal asymptotic behavior quite hard, 
requiring an effective control of the scaling corrections. For this purpose we exploit 
improved Hamiltonians, for which the leading scaling corrections are suppressed for 
any thermodynamic quantity, and improved observables, for which the leading scaling 
corrections are suppressed for any model belonging to the same universality class. 

The results of the finite-size scaling analysis provide strong numerical evidence 
that phase transitions in three-dimensional randomly site-diluted and bond-diluted 
Ising models belong to the same randomly dilute Ising universality class. We obtain 
accurate estimates of the critical exponents, v = 0.683(2), ij = 0.036(1), a — —0.049(6), 
7 = 1.341(4), f3 = 0.354(1), S = 4.792(6), and of the leading and next-to-leading 
correction-to-scaling exponents, uj = 0.33(3) and uj2 — 0.82(8). 
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1. Introduction 

The effect of quenclied random disorder on tlie critical behavior of Ising systems 
has been much investigated experimentally and theoretically. Typical physical 
examples are randomly dilute uniaxial antiferromagnets, for instance, FepZni_pF2 
and MnpZni_pF2, obtained by mixing a uniaxial antiferromagnet with a nonmagnetic 
material. Experiments (see, e.g., Ref. [1] for a review) find that, for sufficiently low 
impurity concentration 1 — p, these systems undergo a second-order phase transition 
at T(.{p) < Tc{p = 1). The critical behavior appears approximately independent of the 
impurity concentration, but definitely different from the one of the pure system. These 
results support the existence of a random Ising universality class which differs from the 
one of pure Ising systemsjj] 

A simple lattice model for dilute Ising systems is provided by the three-dimensional 
randomly site-diluted Ising model (RSIM) with Hamiltonian 

Hs = -J '^PxPyCr^ay, (1) 

<xy> 

where the sum is extended over all nearest-neighbor sites of a simple cubic lattice, are 
Ising spin variables, and px are uncorrelated quenched random variables, which are equal 
to one with probability p (the spin concentration) and zero with probability 1 — p (the 
impurity concentration). Another related lattice model is the randomly bond-diluted 
Ising model (RBIM) with Hamiltonian 

Ti-b = -J ^ jxy CTxay, (2) 

<xy> 

where the bond variables jxy are uncorrelated quenched random variables, which are 
equal to one with probability p and zero with probability 1 — p. Note that the RSIM 
can be seen as a RBIM with bond variables j^y = PxPy But in this case the bond 
variables are correlated. For example, the connected average of the bond variables 
along a plaquette does not vanish as in the case of uncorrelated bond variables: indeed 
(note that j^y = p^) 



\{(jxy-Jxy)=p'-¥ + ^p' -p'. (3) 

□ 

Above the percolation threshold of the spins, these models undergo a phase transition 
between a disordered and a ferromagnetic phase. Its nature has been the object of 
many theoretical studies, see, e.g., Refs. [H El [6l [TJ [8] for reviews. A natural scenario 
is that the critical behavior of the RSIM and the RBIM, at any value of p above the 
spin percolation threshold, belongs to the same universality class, which we will call 
randomly dilute Ising (RDIs) universality class. 

I Unixial antiferromagnets undergo a phase transition also in the presence of a magnetic field H . In 
the absence of dilution, for small H the transition is in the Ising universality class: the magnetic field 
does not change the critical behavior. In the presence of dilution instead, H is relevant and the critical 
behavior for H ^ belongs to the random- field Ising universality class j2j . The crossover occurring for 
— > is studied in Ref. [3|. 
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Figure 1. RG trajectories in the renormalized coupling {u,v) plane starting from 
the Gaussian FP, for several values of the ratio s = uq/vq (from Ref. [5]). We also 
report the RG trajectory connecting the Ising FP with the RDls FP, relevant for the 
Ising-to-RDIs crossover. 



The RDIs universality class can be investigated by field-theoretical (FT) methods, 
starting from the Landau- Ginzburg- Wilson 0^ Hamiltonian 




where 0j is an A^-component field. By using the standard replica trick, it can be 
shown that, for uq < and in the limit — > 0, this model corresponds to a system 
with quenched disorder effectively coupled to the energy density. Fig. [1] shows the 
renormalizat ion-group (RG) flow [9] in the u, v plane, where u, v are the renormalized 
couplings associated with the Hamiltonian parameters uq and vq. The RG flow has a 
stable flxed point (FP) in the region u < 0, which attracts systems with —1 ^ uq/vq < 0. 
The standard Ising FP at m = is unstable, with a crossover exponent [1] = ais, 
where [10] ajs = 0.1096(5) is the specific-heat exponent of the Ising universality 
class. The stable RDIs FP determines the critical behavior at the phase transition 
between the disordered and the ferromagnetic phase that occurs in RDIs systems. 
Therefore, it is expected to determine the critical behavior of the RSIM and of the 
RBIM above the spin percolation threshold. The critical exponents at the RDIs FP 
have been computed perturbatively to six loops in the three-dimensional massive zero- 
momentum scheme |TT1[T2]. Even though the perturbative series are not Borel summable 
[T3l|lll[T5], appropriate resummations provide quite accurate results [11]: z/ = 0.678(10), 
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a = -0.034(30), r] = 0.030(3), 7 = 1.330(17), P = 0.349(5). Moreover, the leading 
scaling corrections are characterized by a small exponent u = 0.25(10), which is much 
smaller than that occurring in pure Ising systems, u ^ 0.8 [Ml E]. Experiments find 
[T| z/ = 0.69(1), a = —0.10(2), and /? = 0.350(9), which are in reasonable agreement 
with the FT estimates (there is only a small discrepancy for a). We also mention that 
approximate expressions for the RDIs critical equation of state have been reported in 



Monte Carlo (MC) simulations of the RSIM and the RBIM have long been 
inconclusive in settling the question of the critical behavior of these models. While 
the measured critical exponents were definitely different from the Ising ones, results 
apparently depended on the spin concentration, in disagreement with RG theory. The 
question was clarified in Ref. [18], where the apparent violations of universality were 
explained by the presence of large concentration-dependent scaling corrections, which 
decay very slowly because of the small value of the exponent uj, uj = 0.37(6). Only if 
they are properly taken into account, the numerical estimates of the critical exponents 
of the RSIM become dilution-independent as expected. Ref. [18] reported the estimates 
u = 0.6837(53) and rj = 0.0374(45), which are in agreement with the FT results. These 
results were later confirmed by MC simulations [19] of the RSIM at p = 0.8, which is 
the value of p where the leading scaling corrections appear suppressed [HI [20]. A FSS 
analysis of the data up to lattice size L = 256 gave z/ = 0.683(3) and 77 = 0.035(2). On 
the other hand, results for the RBIM have been less satisfactory. Recent works, based 
on FSS analyses [21] of MC data up to L = 96 and high-temperature expansions [22], 
have apparently found asymptotic power-law behaviors that are quite dependent on the 
spin concentration. Such results may again be explained by the presence of sizeable and 
spin-concentration dependent scaling corrections. 

The RSIM and the RBIM, and in general systems which are supposed to belong 
to the RDIs universality class, are examples of models in which the presence of slowly- 
decaying scaling corrections makes the determination of the asymptotic critical behavior 
quite difficult. In these cases, the universal critical behavior can be reliably determined 
only if scaling corrections are kept under control in the numerical analyses. For example, 
the Wegner expansion of the magnetic susceptibility x is generally given by [23] 




where t = 1 — /3//?c is the reduced temperature. We have neglected additional terms due 
to other irrelevant operators and terms due to the analytic background present in the 
free energy [2^ [25l [26] . In the case of the three-dimensional RDIs universality class we 
have [l8l[TT] A = uju ~ 0.2, which is very small, and [1^ A2 = a;2i^ ~ 0.5. Analogously, 
the finite-size scaling (FSS) behavior at criticality is given by [25l [27] 



Ref. [IT]. 




X = Ct~^ (1 + ao,it + ao,2t^ + ... + ai,it^ + 01,2^^"^ + . . . 

+ &i,2t'+'^ + . . . + a2,it^^ + ...), 



(5) 



X = cL^-"" (1 + ai,iL-" + ai,2i^^'" + • • • + a2,ii^""' + ■■■), 



(6) 



where ui ^ 0.3 and 002 ~ 0.8 for the RDIs universality class. 
The main purposes of this paper are the following: 
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(i) We wish to improve the numerical estimates of the critical exponents associated 
with the asymptotic behavior and with the leading scaling corrections. 

(ii) We wish to provide robust evidence that the critical behaviors of the RSIM and of 
the RBIM belong to the same RDIs universality class, independently of the impurity 
concentration. 

For these purposes, we perform a high-statistics MC simulation of the RSIM for 
p = 0.8 and p = 0.65, and of the RBIM for p = 0.7 and p = 0.55, for lattice sizes up 
to L = 192. The critical behavior is obtained by a careful FSS analysis, in which a RG 
invariant quantity (we shall use = ^/L) is kept fixed. This method has significant 
advantages [2HI ESI EZ] with respect to more standard approaches, and, in particular, 
it does not require a precise estimate of Pc- Our main results can be summarized as 
follows. 

• We obtain accurate estimates of the critical exponents: u = 0.683(2) and rj = 
0.036(1). Then, using the scaling and hyperscaling relations a = 2— 3z/, 7 = (2— 77)//, 
P = + ri)/2 and 5 = {5 - r/)/(l + t]), we obtain a = -0.049(6), 7 = 1.341(4), 
p = 0.354(1), and 5 = 4.792(6). 

• We obtain accurate estimates of the exponents associated with the leading scaling 
corrections: u = 0.33(3) and UJ2 = 0.82(8). Correspondingly, we have A = = 
0.22(2) and A2 = cu2Z/ = 0.56(5). 

• For both the RSIM and the RBIM we estimate the value p* at which the leading 
scaling corrections associated with the exponent u vanish for all quantities. For the 
RSIM, we find p* = 0.800(5). This result significantly strengthens the evidence that 
the RSIM with p = 0.8 is improved, as already suggested by earlier calculations 
[IHIIIH]. For the RBIM, we find p* = 0.54(2). 

• We provide strong evidence that the transitions in the RSIM and in the RBIM 
belong to the same universality class. For this purpose, the knowledge of the leading 
and next-to-leading scaling correction exponents is essential. We also make use of 
improved observables characterized by the (approximate) absence of the leading 
scaling correction for any system belonging to the RDIs universality class. 

The paper is organized as follows. In Sec. [2] we report the definitions of the 
quantities which are considered in the paper. In Sec. [3] we summarize some basic FSS 
results which are needed for the analysis of the MC data. In Sec. Hlwe give some details 
of the MC simulations. Sec. [5] describes the FSS analyses of the MC data of the RSIM 
ai p = 0.8 and p = 0.65 at fixed R^, which lead to the best estimates of the critical 
exponents. FSS analyses at fixed P of the RSIM at p = 0.8 are presented in Sec. [61 
Finally, in Sec. [7|we analyze the data for the RBIM at p = 0.55 and p = 0.7, and show 
that the RBIM belongs to the same universality class as the RSIM. In Appendix A we 



determine the next-to-leading correction-to-scaling exponent by a reanalysis of the 
FT six-loop perturbative expansions reported in Ref . [30] . In Appendix B we discuss the 
problem of the bias in MC calculations of disorder averages of combinations of thermal 
averages. 
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2. Notations 

We consider Hamiltonians ([T]) and ([2]) with J = 1 on a finite simple-cubic lattice with 
periodic boundary conditions. In the case of the RSIM, given a quantity O depending 
on the spins {a} and on the random variables {p}, we define the thermal average at 
fixed distribution {p} as 

where Z({p}) is the sample partition function. Then, we average over the random 
dilution considering 



{Om = J[dp]{0){(3,{p}), (8) 

where 

[dp]=ll[p6ip,-l) + il-p)5ip,)]. (9) 

Analogous formulas can be written for the RBIM, taking 

m = n iP^ij-y -!) + (!- P)^ij^y)]- (10) 

<xy> 

We define the two-point correlation function 

G{x) = {poao p^a^) (RSIM), 

G{x)=J^^ (RBIM), (11) 

the corresponding susceptibility x. 



X 



^G(x), (12) 



and the correlation length ^, 

^2 ^ g(0) - g(9min) ^ ^^3^ 
^min^(5'min) 

where q^i^i = (Stt/L, 0, 0), g = 2 sing/2, and G(g) is the Fourier transform of G{x). 

We also consider quantities that are invariant under RG transformations in the 
critical limit. We call them phenomenological couplings and generically refer to them 
by using the symbol R. Beside the ratio 

= (14) 

we consider the quartic cumulants and U22 defined by 

f/4 - 5, ^ ^^i^, (15) 

where 

Pfc = ((^p,a, )'=) (RSIM), (16) 

X 

Pk = {{Y,a^f) (RBIM). 
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We also define their difference 

f/rf = f/4 - f/22. (17) 

Finally, we consider the derivative of the phenomenological couplings R with respect to 
the inverse temperature, i.e., 

(18) 

which allows one to determine the critical exponent v. 

3. Finite-size scaling 
3.1. General results 

In this section we summarize some basic results concerning FSS, which will be used 
in the FSS analysis. We consider a generic RDIs model in the presence of a constant 
magnetic field H and in a finite volume of linear size L, and the disorder-averaged 
free-energy density 



H(3,H,L) = j^\nZ{(3,H,L), (19) 

where d is the space dimension (d = 3 in our specific case). In analogy with what is 
found in systems without disorder (see, e.g., Refs. [251 EI]), we expect the free energy 
to be the sum of a regular part jFi,gg(/3, if, L) and of a singular part ^sing(/5, L). The 
regular part is expected to depend on L only through exponentially small terms, while 
the singular part encodes the critical behavior. The scaling behavior of the latter is 
expected to be: 

^smg(Mt, Uh, {ui}, L) = L-''j^,i^^{Ly'ut, ly^Uh, {ly^Ui}), (20) 

where ut = ui, Uh = U2, {ui} with z > 3 are the scaling fields associated respectively 
with the reduced temperature t {ut ~ t), the magnetic field H {uh ~ H), and the other 
irrelevant perturbations with ?/j < 0. They are analytic functions of the Hamiltonian 
parameters — in particular, of t and H — and are expected not to depend on the linear 
size L [25]. Since Ut and Uh are assumed to be the only relevant scaling fields, yi < 
for i > 3. Thus, in the infinite- volume limit and for any t, the arguments L^'Ui go to 
zero. One may thus expand J^sing with respect to L^^Ui obtaining all scaling corrections. 
The RG dimensions of the relevant scaling fields Ut and Uh are related to the standard 
exponents u and r] hj yt = 1/z/ and yh = {d -\- 2 — ri)/2. The correction-to-scaling 
exponents uj and 002 introduced in Sec. [T] are related to the RG dimensions of the two 
leading irrelevant scaling fields: u = —y^ and UJ2 = —Vi- 

The scaling behavior of zero-momentum thermodynamic quantities can be obtained 
by performing appropriate derivatives of Eq. (!20!) with respect to t and H. For instance, 
for the susceptibility at if = we obtain 



xiP,L) = k^L 



2- 



j>3 k 



+ Xreg(/?), (21) 
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where we have neglected terms scahng as L^^^'^y^^ j^yt-Vh^ L^^+vj-'^yh^ which arise 
from the H dependence of Ut and Ui (see, e.g., Refs. [211 [32] for a discussion in the infinite- 
volume limit). The functions Xo{^) ^i-nd Xi,k{z) are smooth, finite for z — 0, and universal 
once one chooses a specific normalization condition (which must be independent of the 
Hamiltonian parameters) for the scaling fields and for the susceptibility (which amounts 
to properly choosing the model-dependent constant k^). The function XregiP) represents 
the contribution of the regular part of the free-energy density and is L independent 
(apart from exponentially small terms). For t ^ we have Ut{t = 0) = 0, while, 
generically, we expect Ui(t = 0) 7^ 0. Expanding Eq. (pT!) for L — > 00, one obtains 
Eq. (ED. 

Analogous formulae hold for the 2n-point susceptibilities. They allow us to derive 

the scaling behavior of the quartic cumulant R = U4. We obtain in the FSS limit 

00 

R{(], L) = roiutLy^) "-iA^tiy') ^'L''' + ^reg(/?), (22) 

i>3 k 

where again several irrelevant terms have been neglected. As before, the functions ro{z) 
and Ti^k^z) are smooth, finite for 2; — 0, and universal. The function rreg(/?) is due to 
the regular part of the free energy and gives rise to scaling corrections of order L''"^. 
For t — >■ 0, writing Ut = Cft + O(t^), we can further expand Eq. fl22l) . obtaining 

R{P, L) = R* + r[,(0)Q tLy^ + rs,o{0)u,{t = 0) L^^ + r^M^^it = Of L^y^ 

+ n,o{o)u4{t = 0) ly* + 0{fL^y\L^y\tLy'+y\L^y\ l^^), (23) 

where R* = ro(0) and we have not written explicitly the corrections due to rreg(/3). Note 
that no analytic 1/L corrections are expected [25l [27] . 

The scaling behavior of U22 can be derived analogously. In this case, one should 
start from the two- replica free-energy density (see, e.g., Ref. |T3|, App. B) 

^2(/5, H,, H2, L) = ^lnZ{(3,H,,L)lnZ{(3,H2,L) (24) 

and assume, as before, that J^2{P, Hi, H2, L) is the sum of a regular part and of a singular 
part. The latter should scale as 

^2,sing(Mt, Mhi, Mh2, {^i}, L) = L-'^J^singiL^Ut, ly^Uh,, ly^Uh,, {L^'Mi}) , (25) 

where we have now two magnetic scaling fields such that Uhi ~ Hi and ~ H2. 
Taking the appropriate derivatives, one can verify that Eq. (12^ holds for U22 too. 

The discussion we have presented applies only to zero-momentum quantities. In 
order to derive the scaling behavior of the correlation length, we should also consider 
quantities that are defined in terms of the field variables at nonzero momentum. They 
can be derived from a free-energy density associated with a model in which the magnetic 
field is site dependent. The general analysis by Wegner [23] indicates that Eq. f[25l) 
should still hold (at least in systems without disorder; we assume here that the same 
holds for quenched averages): the presence of a site-dependent H only modifies the 
scaling fields that become functionals of H. For these reasons, we conjecture that also 
^ has an expansion similar to ( l22l) : in particular, we expect corrections proportional 



Universality class of 3D site- diluted and bond-diluted Ising systems 



9 



to u^L~^'^ and u'lL~'"^^. The momentum dependence of the scahng fields will give rise 
to additional terms, the leading one being proportional to (for a discussion in the 
two-dimensional Ising model, see Ref. [33j, p. 8161). Additional 1/L^ corrections arise 
from our particular definition of ^ [33] . The analyses we shall report below confirm this 
conjecture. 

The scaling fields Ui depend on the model. Thus, if one considers families of models 
that depend on an irrelevant parameter, by a proper tuning one can find Hamiltonians 
for which us{t = 0) = 0. In this case, in the FSS limit t — 0, L ^ oo at fixed tL^*, 
all corrections proportional to L*^^^ = L^'^'^ vanish. Note that, since the scaling field 
depends only on the model, such a cancellation occurs for any quantity one considers. 
In all cases the leading correction-to-scaling term behaves as L~'^^. Models such that 
W3(t = 0) = will be called improved models. Since the leading scaling correction 
vanishes, one should observe a faster approach to the scaling limit. 

In this paper we shall also consider improved observables. An improved 
phenomenological coupling is such that ^3,0(0) = 0. As a consequence, in the FSS limit 
at fixed UfLy* = it does not show leading scaling corrections proportional to L^*^. Note 
that this cancellation is only observed on the line in the (t, L) plane such that tL^* = 0, 
but not in the generic FSS limit. As a consequence, if R is improved, its derivative R' 
is not. Note also that, while in improved models all corrections proportional to L~^^ 
vanish, here only the leading one vanishes: corrections proportional to L~^^ are still 
present. Improved observables satisfy a very important property. Since the function 
r3fi{z) is universal, the cancellation occurs in any model|§| 

The thermal RG exponent yt = 1/z/ is usually computed from the FSS of the 
derivative R' of a phenomenological coupling R with respect to /3 at Pc- Using Eq. (!22ll 
one obtains 



1 + asusit = 0)Ly-' + a^Uiit = 0)1^^ + h^it = 0) 1^'"^' 



(26) 



where Sq, 03, 04 are constants. The leading correction scales as L^'^ = L^'^. In improved 
models, in which us{t = 0) = 0, the leading correction is of order L^* = L""^^. Note 
that corrections proportional to Ly^~y* = j^-^-'^/'^ are still present even if the model is 
improved. 

§ In principle one could define an improved observable along any line with fixed tL^* . Such observables 
would be improved in any model on the same line utL^* . However, these observables are not very useful 
in practice. Suppose indeed that one determines an improved observable in a model (call it model A) 
along the line such that tL^* — kA- li ut = CAt, this means that improvement is observed on the line 
UtLy* = UaCa- Then consider a second model (model B). In model B the observable is improved on 
the same line = UaCa- Now ut = est with cb 7^ ca- Thus, in model B improvement is observed 

on the line iL^' = kACA/cs- Thus, in order to use the improved observable in model B one should 
determine the model-dependent ratio ca/cb, which is quite inconvenient. No such difficulty arises for 
fcyi = 0. All these problems are avoided by using FSS at fixed phenomenological couplings, see next 
subsection. 
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3.2. Finite-size scaling at a fixed phenomenological coupling 

Instead of computing the various quantities at fixed Hamiltonian parameters, one may 
study FSS keeping a phenomenological coupling R fixed at a given value Rf [28]. This 
means that, for each L, one considers (3f{L) such that 

Ri^ = ^fiL),L) = Rf. (27) 

All interesting thermodynamic quantities are then computed at /5 = (3f{L). The 
pseudocritical temperature (3f{L) converges to /5c as L — * oo. The value Rf can be 
specified at will, as long as Rf is taken between the high- and low-temperature fixed- 
point values of R. The choice Rf = R* (where R* is the critical-point value) improves 
the convergence of /5/ to /5c for L ^ oo; indeed [281 EH] Pf — Pc = 0{L~^/'') for generic 
values of i?/, while Pf — Pc = 0{L^^^^~^) for Rf = R*. This method has several 
advantages. First, no precise knowledge of Pc is needed. Secondly, for some observables, 
the statistical error at fixed Rf is smaller than that at fixed P = Pe- 
lf Pf is determined as in Eq. (1271) . one may consider the value of other 
phenomenological couplings Ra at Pf, defining 

Ro^{L) = R^[Pf{L),L]. (28) 

The large-L limit of Ra is universal but depends on Rf (it differs from the critical 
value i?*, unless Rf = R*). Indeed, neglecting scaling corrections in Eq. ( l22l) . we have 
Ra = fafliutW^) and R = r^iutL'^^), where Ta-^Q^z) and tq^z) are universal functions. 
Fixing R = Rf corresponds to fixing a particular trajectory in the t, L plane given by 
UtL^* = Zf, where Zf is the solution of the equation Rf = rQ^Zf). Along this trajectory 
Ra ~ Ra = ra;o{zf), which shows the universality of Ra- 

We can define improved observables also when considering FSS at fixed Rf. Such 
observables show a faster convergence since the corrections to the scaling limit scale 
as L~'^^ and L"^"^. Moreover, at variance with the case discussed in Sec. 13.11 where it 
was practical to define an improved observable only on the line tL^* = 0, here one can 
choose any value for Rf. If one determines an improved observable in a given model for 
a chosen value Rf, this observable is improved in any other model at the same value 
Rf. 

4. Monte Ccirlo simulations 

We performed MC simulations of the RSIM at p = 0.8 and p = 0.65 and of the RBIM 
at p = 0.55 and p = 0.7 close to the critical temperature on cubic lattices of size L^, 
with L = 8, 12, 16, 24, 32, 48, 64, 96, 128. For the RSIM at p = 0.8 we also performed 
simulations with L = 192. For each lattice size, we collected Ns disorder samples, Ns 
varying between 10^ and 6.4 x 10^. In Tables (H [21 [3l HI and [5] we report some details of 
the MC simulations and the results at fixed R^ = ^/L = 0.5943. 

In the simulations we used a combination of Metropolis, Swendsen-Wang (SW) 
cluster [36], and Wolff single-cluster [37] updates, to achieve an effective thermalization 
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Table 1. Run parameters and estimates of /?/ and x fixed — 0.5943 for the 
RSIM at p = 0.8. Ns is the number of disorder samples divided by 1000. 



L 


Ns 






Pf 


X 


8 


200 





285744 


0.286020(24) 


63.391(9) 


12 


200 





285744 


0.285851(14) 


141.581(19) 


16 


200 





285744 


0.285800(6) 


249.910(25) 


24 


100 





285744 


0.285765(6) 


555.67(7) 


32 


100 





285744 


0.285751(4) 


979.36(13) 




100 





285761 


0.285753(4) 


979.27(12) 


48 


60 





285744 


0.2857515(24) 


2173.7(3) 




106 





285748 


0.2857459(18) 


2173.8(3) 




60 





285751 


0.285743(3) 


2174.0(3) 


64 


60 





285742 


0.2857478(18) 


3827.1(6) 




63 





285744 


0.2857443(17) 


3827.4(6) 


96 


30 





285744 


0.2857441(13) 


8491.6(1.9) 


128 


20 





285743 


0.2857428(9) 


14950(4) 




20 





285744 


0.2857438(10) 


14945(4) 


192 


10 





285743 


0.2857430(7) 


33150(12) 




10 





285744 


0.2857435(8) 


33141(13) 



of short- and long-range modes. For each disorder sample, we started from a random 
spin configuration, then we typically performed 300 thermalization steps, each step 
consisting in 1 SW update, 1 Metropolis update, and L single-cluster updates. Then, 
we typically performed 400 measures for lattice sizes L < 64 and 600 measures for larger 
lattices. Between two measurements we usually performed 1 SW update and 2L single- 
cluster Wolff updates. We did several tests of thermalization, performing some runs 
with a larger number of thermalization steps; to test the independence of the results on 
Njn, we also did few runs with a larger number of measures per disorder configuration. 

In the determination of the averages over disorder one should take care of the 
bias that occurs because of the finite number of measures at fixed disorder [38]. A 
bias correction should be introduced whenever one considers the disorder average of 
combinations of thermal averages and, in particular, whenever a reweighting of the 



data is performed. Details are reported in [Appendix B[ Errors are computed from the 



sample-to-sample fluctuations and are determined by using the jackknife method. 

The MC simulations that we present here took approximately 10 CPU years of a 
workstation equipped with an AMD Opteron Processor 246 (2 GHz clock frequency). 

5. Finite-size scaling analysis at fixed of the randomly site-diluted Ising 
model 



In this section we describe the FSS analyses of the MC simulations for the RSIM at 
p = 0.8 and p = 0.65. The analyses are performed at a fixed value of R^. 
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Table 2. MC results for the phenomenological couplings at fixed = 0.5943 for the 
RSIM atp = 0.8. 



L 




U22 


U4. 


Ud 




R'e 


-U'4 


8 


0.285744 


0.1504(6) 


1.6066(5) 


1.45624(19) 


1.8021(14) 


18.250(9) 


32.153(20) 


12 


0.285744 


0.1490(4) 


1.6197(3) 


1.47067(13) 


1.8134(9) 


32.861(14) 


58.99(3) 


16 


0.285744 


0.1484(5) 


1.6262(4) 


1.47779(14) 


1.8192(10) 


49.985(22) 


90.40(5) 


24 


0.285744 


0.1480(6) 


1.6328(5) 


1.48474(21) 


1.8252(14) 


90.24(6) 


164.51(14) 


32 


0.285744 


0.1482(6) 


1.6362(5) 


1.48798(19) 


1.8289(13) 


137.28(8) 


251.20(20) 




0.285761 


0.1495(7) 


1.6370(6) 


1.48752(22) 


1.8313(15) 


137.12(9) 


251.23(19) 


48 


0.285744 


0.1476(8) 


1.6394(7) 


1.49179(25) 


1.8313(17) 


248.29(19) 


456.2(4) 




0.285748 


0.1494(6) 


1.6407(5) 


1.49132(20) 


1.8348(13) 


248.20(14) 


456.3(3) 




0.285751 


0.1497(9) 


1.6409(8) 


1.49117(27) 


1.8356(20) 


248.04(20) 


456.4(5) 


64 


0.285742 


0.1472(9) 


1.6410(7) 


1.49380(26) 


1.8324(18) 


378.5(3) 


695.8(8) 




0.285744 


0.1475(8) 


1.6412(7) 


1.49365(26) 


1.8330(17) 


378.7(3) 


696.2(7) 


96 


0.285744 


0.1478(12) 


1.6429(10) 


1.4950(4) 


1.8351(25) 


684.2(8) 


1259.9(1.9) 


128 


0.285743 


0.1477(15) 


1.6440(13) 


1.4963(5) 


1.836(3) 


1043.2(1.5) 


1923(4) 




0.285744 


0.1471(18) 


1.6434(15) 


1.4963(5) 


1.835(4) 


1042.7(1.6) 


1921(4) 


192 


0.285743 


0.1486(19) 


1.6459(16) 


1.4973(6) 


1.839(4) 


1889(4) 


3493(12) 




0.285744 


0, 1490(22) 


l.(i4()0(18j 


1.4970(7] 


1. 840(5) 


1887(4] 


3485(14) 



Table 3. MC results at fixed R^ = 0.5943 for the RSIM at p = 0.65. Ns is the 
number of samples divided by 1000. 



L 




^run 


/5/ 


X 


U22 


U4 




R'^ 


8 


200 


0.373250 


0.37372(6) 


50.964(9) 


0.2113(7) 


1.5609(6) 


1.8356(14) 


9.882(13) 


12 


200 


0.371650 


0.37180(3) 


113.804(16) 


0.2016(6) 


1.5751(5) 


1.8372(12) 


17.581(21) 


16 


200 


0.371050 


0.371132(19) 


200.92(3) 


0.1956(6) 


1.5827(5) 


1.8369(13) 


26.371(15) 




200 


0.371347 


0.371080(19) 


201.045(23) 


0.1959(6) 


1.5831(5) 


1.8378(13) 


26.360(15) 


24 


50 


0.370500 


0.370626(21) 


446.94(10) 


0.1896(11) 


1.5926(10) 


1.8390(24) 


46.69(5) 




50 


0.370600 


0.370606(19) 


446.87(13) 


0.1892(12) 


1.5929(10) 


1.839(3) 


46.78(5) 




50 


0.370700 


0.370620(17) 


446.98(9) 


0.1904(11) 


1.5934(10) 


1.8410(24) 


46.66(6) 


32 


100 


0.370420 


0.370457(9) 


787.67(13) 


0.1856(8) 


1.5984(7) 


1.8397(17) 


71.13(6) 




100 


0.370482 


0.370420(9) 


788.01(12) 


0.1853(9) 


1.5982(7) 


1.8392(18) 


71.11(6) 


48 


30 


0.370290 


0.370310(8) 


1749.2(4) 


0.1821(14) 


1.6068(11) 


1.844(3) 


124.59(17) 




30 


0.370304 


0.370304(8) 


1750.0(5) 


0.1823(15) 


1.6067(13) 


1.844(3) 


124.39(15) 




30 


0.370329 


0.370304(9) 


1751.2(5) 


0.1822(15) 


1.6061(12) 


1.843(3) 


124.37(15) 


64 


30 


0.370156 


0.370249(5) 


3079.8(9) 


0.1778(15) 


1.6096(13) 


1.841(3) 


188.1(6) 




30 


0.370249 


0.370243(6) 


3084.6(8) 


0.1761(15) 


1.6081(12) 


1.837(3) 


187.6(6) 


96 


20 


0.370156 


0.370209(4) 


6839.3(2.1) 


0.1742(15) 


1.6149(13) 


1.841(3) 


335.8(1.0) 




20 


0.370209 


0.370215(3) 


6842.2(2.3) 


0.1734(19) 


1.6142(15) 


1.840(4) 


334.4(1.2) 


128 


10 


0.370196 


0.370196(3) 


12059(5) 


0.170(3) 


1.6156(22) 


1.837(6) 


504.0(2.2) 




10 


0.370209 


0.370196(3) 


12059(5) 


0.1702(25) 


1.6162(21) 


1.837(5) 


504.8(2.0) 
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Table 4. MC results at fixed = 0.5943 for tlie RBIM at p = 0.55. Ns is tlie 
number of samples divided by 1000. 



L 


Ns 


Prun 


Pf 


X 


U22 


U4 






8 


640 


0.432500 


0.432539(17) 


79.781(4) 


0.12857(22) 


1.60385(18) 


1.7710(5) 


11.270(3) 


12 


640 


0.432340 


0.432384(10) 


178.251(9) 


0.13629(23) 


1.61874(19) 


1.7959(5) 


20.143(5) 


16 


640 


0.432330 


0.432335(6) 


314.646(17) 


0.1393(3) 


1.62592(20) 


1.8070(5) 


30.544(7) 


24 


256 


0.432340 


0.432303(5) 


699.72(6) 


0.1432(4) 


1.6338(3) 


1.8199(8) 


55.036(20) 


32 


100 


0.432300 


0.432295(5) 


1232.80(18) 


0.1440(7) 


1.6371(6) 


1.8243(14) 


83.76(5) 


48 


62 


0.432293 


0.432294(5) 


2737.8(7) 


0.1443(12) 


1.6401(10) 


1.828(3) 


151.64(16) 


64 


30 


0.432285 


0.432290(4) 


4820.2(9) 


0.1443(10) 


1.6411(8) 


1.8288(21) 


230.27(20) 


96 


10 


0.432279 


0.432291(4) 


10693(3) 


0.1446(20) 


1.6429(16) 


1.831(4) 


417.9(7) 




on 
ZU 


0.432288 


0.432290(3) 


10693.6(2.4) 


0.1480(14) 


1.6455(11) 


1.838(3) 


416.6(5) 




on 


0.432294 


0.4322901(17) 


18823(4) 


0.1472(14) 


1.6457(12) 


1.837(3) 


635.4(7) 






Table 5. MC results at fixed = ( 


3.5943 for the RBIM at p = 0.7 


. iV, is the number 






of samples divided by 1000. 














Pvun 


Pf 


X 


U22 


U4 


C/im 




8 


640 


0.325900 


0.325905(9) 


79.519(3) 


0.08327(15) 


1.61733(12) 


1.7256(3) 


18.308(3) 


12 


640 


0.326200 


0.326253(5) 


177.657(8) 


0.09027(17) 


1.63479(13) 


1.7521(3) 


33.299(6) 


16 


640 


0.326400 


0.326409(3) 


313.629(13) 


0.09444(17) 


1.64355(14) 


1.7663(4) 


51.082(8) 


24 


250 


0.326560 


0.326547(3) 


697.55(5) 


0.0991(3) 


1.65209(24) 


1.7810(6) 


93.642(24) 


32 


100 


0.326600 


0.326605(3) 


1229.06(14) 


0.1021(5) 


1.6564(4) 


1.7891(10) 


144.08(6) 


48 


30 


0.326660 


0.326651(3) 


2729.0(5) 


0.1059(9) 


1.6606(7) 


1.7983(18) 


264.59(23) 


64 


30 


0.326629 


0.3266746(16) 


4802.1(8) 


0.1073(9) 


1.6618(8) 


1.8013(20) 


408.1(8) 




30 


0.326640 


0.3266792(19) 


4802.9(7) 


0.1069(8) 


1.6618(7) 


1.8008(17) 


408.2(6) 


96 


20 


0.326649 


0.3266919(12) 


10650.6(1.9) 


0.1117(11) 


1.6645(10) 


1.8097(24) 


749.4(2.6) 




20 


0.326663 


0.3266922(13) 


10650.4(2.2) 


0.1122(12) 


1.6655(10) 


1.811(3) 


749.0(2.2) 


128 


10 


0.326698 


0.3266982(11) 


18756(5) 


0.1144(18) 


1.6665(15) 


1.815(4) 


1151(3) 




10 


0.326706 


0.3266977(12) 


18775(5) 


0.1131(17) 


1.6646(14) 


1.812(4) 


1150(3) 



5.1. FSS at fixed phenomenological coupling R 

Instead of performing the FSS analysis at fixed Hamiltonian parameters, we analyze 
the data at a fixed value i?/ of a given phenomenological coupling R, as discussed in 
Sec. 13. 2[ The most convenient choice for the value Rf is Rf ~ R*, where R* is the 
asymptotic value of R at /5c- Data at fixed R = Rf are obtained by computing i?(/3) in 
a neighborhood of (3c- This is done by reweighting the MC data obtained in a simulation 
at P = /5run ~ Pc- Given R{P), one determines the value Pf such that R{P = /?/) = Rf. 
All interesting observables are then measured at their errors at fixed R = Rf are 
determined by a standard jackknife analysis. 

In the following we present the FSS analysis at fixed R^ = ^/L, choosing R^j = 
0.5943, which is the estimate of obtained in Ref. [I9]: i?| = 0.5943(9). In Tables [II [2l 
and [3] we report the results obtained for the RSIM at p = 0.8 and p = 0.65, respectively. 
We also performed FSS analyses at fixed f/4 = U^j = 1.650 [from Ref. [19] that quotes 
= 1.650(9)]. We do not report the corresponding results because they are consistent 
with, though slightly less precise than, those obtained at fixed R^ = 0.5943. 
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Figure 2. MC results tJ22{L) and estimates of C/'|2(^min) as obtained in different fits 
at fixed = 0.5943 for the RSIM at p = 0.8. Some data are slightly shifted along the 
X-axis to make them visible. On the x axis we report L when plotting the MC data 
and the minimum lattice size imin used in the fit when plotting the fit results. The 
dotted lines correspond to the final estimate U22 — 0.148(1). 



FSS at fixed R has the advantage that it does not require a precise knowledge 
of the critical value (3c- But there is another nice side effect: for some observables the 
statistical errors at fixed Rf are smaller than those at fixed (3 (close to (3^. For example, 
in the case of the RSIM at p = 0.8, we find 



err[x|/ 



err[x|i?5=o.5943] 

err[?74|/;J 
err[f/4|i{^ =0.5943] 

err[x'lfe] 
err [x' I ij^ =0.5943] 



10, 



1.7, 



2.0, 



err[x|/ 



err[x|[/4=i.65o] 
err[[/22|/3j 

err[f/22|_Rj=0.5943] 

err[i?^|/3j 



1.7, 



(29) 



err[i?' 



1.1, 



1.^ 



=0.5943 



which are approximately independent of L. Similar numbers are found for the other 
models that we have simulated. 



5.2. Universal values ofU22 and U^ at fixed R^ = 0.5943 

Here we determine the universal large-L limits of the quartic cumulants U22 and f/4 at 
fixed R^ = 0.5943 — we call them U22 and U4, respectively — using the data at p = 0.8 
reported in Table[21 The data of f722 in the range 8 < L < 192 shown in Fig.[2]have a very 
small L dependence. In particular, there is no evidence of scaling corrections associated 
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Figure 3. MC results U4{L) and estimates of C/|(i„iin) as obtained in different fits 
at fixed = 0.5943 for the RSIM at p = 0.8. In the second fit (a + bL'^) e is a free 
parameter. As in Fig.[2l the x axis corresponds to L (MC data) and imin (fit results). 
The dotted lines correspond to the final estimate — 1.648(3). 



with the leading correction-to-scahng exponent lo ^ 0.3. This confirms ear her results 
p!8l |20[ [19] indicating that leading scahng corrections in the RSIM at p = 0.8 are very 
small. Moreover, also the next-to-leading scaling corrections associated with uj2 ~ 0.8 
are quite small. Indeed, if we fit the data with L > Lmin to a constant, a fit with 
X^/DOF < 1 (DOF is the number of degrees of freedom of the fit) is already obtained 
for Lmin = 12. The corresponding estimates of are independent of L^nin within error 
bars, see Fig. [21 We also fit the data to a -|- bL~'' with e = 0.33 and e = 0.82, which are 
the best MC estimate of u (see Sec. 15.31) and the FT estimate of U2 (see [Appendix A[ ), 
respectively. In both cases b ^ within errors already for Lmin = 16. Finally, we also do 
fits keeping e as a free parameter. The estimates of U22 are independent of Lmin within 
error bars; for Lmin = 8 we get U22 = 0.1483(3), with x^/DOF = 0.8. The estimates 
of U22 obtained in the fits with e = 0.82 are also plotted in Fig. [2] versus the minimum 
lattice size L^^^^ of the data considered in the fits. Collecting results, we obtain the final 
estimate 

U;2 = 0.148(1). (30) 

We perform a similar analysis of f/4. The data are shown in Fig. [31 In this case, there 
is clear evidence of scaling corrections. A fit of the data to a + foL"*^, taking e as a free 
parameter, gives = 1.6472(7) and e = 0.95(5) for Lmm = 8: There is no evidence 
of scaling corrections with exponent u ~ 0.3. This is confirmed by fits to a -|- bL~'^ 
with e = 0.3 fixed. The leading term a varies significantly with L^ain, exactly as the 
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Figure 4. MC results Ud{L) and estimates of f/^(Lmin) as obtained in different fits 
at fixed = 0.5943 for the RSIM at p = 0.8. In the second fit (a + bL'") e is a free 
parameter. As in Fig. [21 the x axis corresponds to L (MC data) and Lmin (fit results). 
The dotted lines correspond to the final estimate — 1.500(1). 



original data (see Fig. [3]). The numerical results for 1/4^ are much better described by 
assuming that the leading scaling corrections are proportional to L^^'^ ~ L^°'^^. A fit to 
a + bL~'^-^^ gives estimates of that show a very tiny dependence on L^in, see Fig. |3l 
The dependence on uj2 is also small: if U2 varies between 0.74 and 0.90 [this corresponds 
to considering the FT estimate uj2 = 0.82(8)], estimates change by less than 0.001 for 
Lmin < 32. These analyses lead to the final estimate 

f/; = 1.648(3). (31) 

We finally consider the difference IJd = U^ — U22, whose data are very precise because 
of an unexpected cancellation of the statistical fiuctuations, see Table [2l In Fig. H] we 
show the results of the same fits as done for IJ4,. They lead to the estimate 

U2 = 1.500(1), (32) 

which is perfectly consistent with the estimates obtained for U22 and tJl- 

5.3. Estimate of the leading correction-to-scaling exponent u 

In order to estimate u, we use the data at p = 0.65. We consider the differences 

A22,a = U22ip = 0.65; L) - U22ip = 0.8; L), (33) 
A,,, = Udip = 0.65; L) - f/,(p = 0.8; L), (34) 
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Figure 5. Estimates of cu from fits of A22,a = U22{p — 0.65) — U22{p = 0.8), A22,b = 
U22{p = 0.65) - Ui2, Ad,a = Ud{p = 0.65) -Ud{p = 0.8), A^,,, = Ud{p = 0.65) - U^, as 
a function of Lmin, which is the minimum lattice size considered in the fits. A22,a and 
A22,6 are fitted as In A = a — uj In L; Ad^a and Ad^b are fitted as In A = a — ujlnL + bL^'^. 
Some data are slightly shifted along the x-axis to make them visible. The dotted lines 
correspond to the final estimate uj = 0.33(3). 



and also 

A22,b = U22{p = 0.65; L) - U;^, (35) 
Ad,b = Udip = 0.65; L) ~ Ul (36) 

where U22 and f/J are given in Eqs. (130|) and (1321) . Because of the universahty of the 
large-L hmit of the phenomenological couphngs — hence, they do not depend on p — they 
behave as 

A ^ CA.llL"" + CA,12i^"'" + ■ ■ ■ + CA,2li^""' + • • • (37) 

The quantities defined in Eqs. (1551) and (15^ are well fitted to CA,ii-Z^~'^[Ji] For instance, 
the analysis of A22,a gives u = 0.342(10) and u = 0.352(24) for L^i^ = 8 and 24, 
respectively; in both cases x^/DOF ^ 0.5. Such a fit instead gives a relatively large 
X^/DOF (x^/DOF = 2.3 for L^i^ = 12) when applied to Ud, essentially because the 
data of Ud have a better relative precision. Moreover, a clear systematic drift is observed 
when varying Lmin- Therefore, we must include the next-to-leading corrections. We thus 
performed fits of the form ca,ii-Z^~'^(1 + dL~'^), where e is an effective exponents that 
takes into account two next-to-leading corrections and should vary in [cj, uj2 — oj\. Given 

II We will often say that we fit a quantity O to aL^ or to aL^{l + 6i^), y <i). What we really do is a 
fit of In O to In a + X In L and to In a + a; In L + bL^ . 
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Figure 6. Plot of U22 versus L'"^ with lj = 0.33 for tlie RSIM (s) at p = 0.8 and 
p = 0.65 and for the RBIM (b) at p 0.7 and p = 0.55. 

the results obtained from the analysis of U22 and the FT estimate uj2 = 0.82(8), we 
have taken e G [0.3,0.6]. The dependence on e is small: for instance, for Lmm = 16, 
the analysis of Ad,a gives u = 0.35(3), 0.33(2), 0.32(2) for e = 0.3, 0.5, 0.6. The results 
corresponding to e = 0.33 and e = 0.49 ^ U2 — uj are reported in Fig. [5] as a function of 
-^min, the smallest lattice size used in the analysis. They lead to the estimate 

uj = 0.33(3), (38) 

which is in agreement with the FT six-loop result [11] u = 0.25(10) (we also mention 
the five-loop result u = 0.32(6) of Ref. [I2]) and with the MC result [H] a; = 0.37(6). 

In writing Eq. (|37|) we assumed that the scaling limit does not depend on p. 
We can now perform a consistenty check, verifying whether U22 and for p = 0.65 
converge to the estimates ( l30l) and ( l3Ti) . A fit of U22 to U22 + £22,11-^"'^ + £22,2-^^"^ gives 
0'*2 = 0.154(4), 0.152(2), 0.149(6), 0.148(5) for (L^^in, e) = (8,2cj), (8,^2), (12,2cj), 
(12,ci;2), respectively. Here we used u = 0.33 as determined above and 002 = 0.82. 
These results are in agreement with the estimate fl30|l . Such an agreement is also clear 
from Fig. [6l where we plot U22 versus L~'^. The same analysis applied to 1/4^ gives 
Ul = 1.640(4), 1.644(3), 1.640(5), 1.644(4), for the same values of {L^in,e), which are 
compatible with the estimate (!3T|) . 

5.4- Determination of the improved RSIM 

We now estimate the value p* of the spin concentration that corresponds to an improved 
model: for p = p* the leading scaling corrections with exponent ut vanish. We already 
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Figure 7. Results of the fits to estimate the amplitude of the leading scaling correction 
in U22- 



know that the RSIM with p = 0.8 is approximately improved, so that p* ~ 0.8. In the 
following we make this statement more precise. We consider U22, which has small next- 
to-leading scaling corrections, and determine the value of p at which the leading scaling 
corrections to this quantity vanish. As we remarked in Sec. 13.11 the same cancellation 
occurs in any other quantity. 

To determine p*, we fit the data of U22 at p = 0.8 andp = 0.65 to tJ 22 + C22,ii{p)L~'' 1 
and the difference A22,a defined in Eq. ( I33l) to ca,ii-^v~'^. We have performed fits in which 
e is fixed, e = uj = 0.33(3), and fits in which it is taken as a free parameter. The results 
of the fits with e = uo are shown in Fig. [7] for several values of L-^ 
size included in the fit. We obtain 



the smallest lattice 



CA,ii = 0.122(6), 

C22,ii(p = 0.65) = 0.122(6), 

C22,ii(p = 0.8) = 0.000(4), 

where the first two estimates correspond to L„ 



16 and the last one to Lr, 



(39) 



48; 



errors are such to include the results of all fits. These results give us the upper bound 



C22,ll(P = 0.8) 



C22,ll(P = 0.65) 



< 



122 



1 

30' 



(40) 



For p = 0.8 scaling corrections are at least a factor of 30 smaller than those occurring 
for p = 0.65. This bound will be useful in the following to assess the relevance of the 
"systematic error" in the fits of the data at p = 0.8 due to possible residual leading 
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scaling corrections. Indeed, the ratio that appears in Eq. (HP]) does not depend on the 
observable. In the notations of Sec. [3], C22,ii(p) = r3 o(-2/)u3(t = 0,p) where r3 o(-2/) is a 
model-independent constant that depends on the observable (in this case on U22)- The 
constant r^Q^zj) drops out in the ratio, since 

C22,ll(pi) ^ u^{t = 0,Pi) 

C22,ll(P2) U3{t = 0,P2)' 

Therefore, given a generic observable 0{p) that behaves as 

0{p) = a(p)L-(l + co,u{p)L-'' + •••), (42) 
we have in all cases 

\co,ii{p = 0.8)/coMp = 0-65)1 < 1/30. (43) 

Estimates (139|) allow us to estimate p*. We obtain p* ^ 0.80. Since, C22,ii(p) = a{p — p*) 
close to p ~ 0.80, the error on p* is simply err[c22,ii(p = 0.8)]/|a|. The constant a — we 
only need a rough estimate since it is only relevant for the error on p* — is determined 

as 

dc22,n 



^ C22,n(p = 0.8) -C22,n(l> = 0.65) ^ 
dp p^p, 0.8-0.65 ^ ^' ^ ^ 

where the reported error is obtained from those of C22,ii(p = 0.8) and C22,ii {p = 0.65). 
This gives 

p* = 0.800(5). (45) 

We are not able to assess the error on a due to the approximation (jS]). Note, however, 
that the dependence on a is small. If we vary a by a factor of 2, p* changes at most by 
0.01. 



5. 5. Determination of the critical temperatures 

We determine the critical temperature by extrapolating the estimates oi (3f reported in 
Tables [T] and [31 According to the discussion reported in Sec. [31 since we have chosen 
= 0.5943 ^ we expect in general that (3f- /3c = ©(L^^/^-^). For p = 0.8, since 
the model is improved, the leading scaling corrections are related to the next-to-leading 
exponent uj2- Thus, in this case Pf — Pc = 0{L~^^^^'^'^) 

In Fig. [31 we show the data for p = 0.8 versus /.^i/'^-'^a taking z/ = 0.68 and 
ui2 = 0.82. The expected linear behavior is clearly observed. A fit to /?c + cL~^^'^~'^^ 
gives = 0.2857429(4). We have also taken into account the error on 002 and on u (the 
error due to the variation of u is essentially negligible compared to the first one). This 
result improves the estimate of (3c obtained in Ref. p[9], i.e. Pc = 0.285744(2). 

For p = 0.65 we fit /?/ to pc + cL-^/^-^ + dL-^/"-' with u = 0.68, cu = 0.33(3), 
e G [0.6,0.9]. We obtain P^ = 0.370174(3). This is consistent with the estimate 
pc = 0.370166(6) {Lmin = 16) reported in Ref. p]. 
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Figure 8. Estimates of /?/ for p = 0.8 versus L-(i/!^+'^2)^ foj. = 0.82 and v = 0.68. 
The dashed Une corresponds to a hnear fit of the data. 



5. 6. Improved phenomenological couplings 

Beside considering improved Hamiltonians — in these particular models any thermody- 
namic quantity does not have leading scaling corrections — one may also consider im- 
proved observables which are such that the leading scaling correction vanishes for any 
Hamiltonian. Here we determine an improved phenomenological coupling by taking an 
appropriate combination of the cumulants and U22, i-e- we consider 

Ui^ = U4 + Cirnf/22. (46) 



In the scaling limit we have generically f/# ~ f/1 + c#,iiL . The constant o 



IS 



determined by requiring Cim^n 

C4,ii(p) 



0, which gives 



(47) 



C22,ll(p) 

where we have written explicitly the p dependence of the coefficients 04^11 (p) and 022,11 (p)- 
As discussed in Sec. 13.1^ the ratio fllTl) is universal within the RDIs universality class 
and, in particular, independent of p. Indeed, 04,11 = ''";74,3,o(^/)w3, 022,11 = ''^c/22,3,o(^/)'W3) 
where r;74^3^o(^) and ru^^^^^^z) are the universal scaling functions defined in Sec. 13.11 
The model-dependent scaling field M3 cancels out in the ratio, proving its universality. 
The combination with the choice ( l47l) has no leading scaling corrections associated 
with the exponent cu, so that, see Sec. EUl t/im = U*^ + ©(L^^"", L-"^''). 
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The ratio be estimated by determining the leading scahng correction 

amphtudes of U22 and [74. Alternatively, one may estimate it from the ratio 

f/22(p = 0.65; L)-f/22(p = 0.8; L) '"^ ^ ' 

Both analyses give consistent results, leading to the estimate 

Cim = 1.3(1). (49) 

Therefore, 

f/im = f/4 + 1.3f/22 (50) 

has (approximately) vanishing leading scaling corrections for any model and any p. 
More precisely, we have the upper bound |cim,ii| ^ 0.1|c22,ii|. Since 04,11 = — CimC22,ii 
andcrf^ii = -(1 + Cim)c22,ii, we also have |cim,ii| < 0.1|c4,ii| and |cim,ii| < 0.05|cd,ii|. The 
leading scaling corrections in tJi^ are at least a factor of 10 smaller than those occurring 
in IJ22 and f/4 and a factor of 20 smaller than those occurring in Ud- This is confirmed by a 
direct analysis of the data of Aim,;,, defined as in Eq. fl35|) . at p = 0.65. A fit to Cim^nL^'^, 
uj = 0.33, gives |cim,ii| ^ 0.005, to be compared with 022,11 = 0.122(6) obtained in 
Sec. [531 In Fig. [9] we show results of fits of tJi^ to U*^ + cL~^ with e = 0.66 ^ 2uj 
and e = 0.82 ^ ^^2. They are consistent with the estimate U-*^^^ = 1.840(4), which can 
be obtained from the estimates of U4, U22, and Ud obtained in Sec. 15. 2[ The improved 
quantity Uirn is particularly useful to check universality, because it is less affected by 
scaling corrections. 

We should note that f/jm is useful in generic models in which corrections are 
generically present, but is not the optimal quantity in improved models in which the 
leading scaling corrections are proportional to C2iL^'^'^ ■ The coefficients C21 can be 
estimated from the data at p = 0.8, obtaining 

C22,21 = 0.01(2), 
C4,21 = -0.21(5), 
Cd,2i = -0.22(5), 

Cim,2i = -0.20(5). (51) 

Since the ratios Ca,2i/cb,2i = ''"a,4,o(2;/)/Tfe,4,o(^/) are universal, the coefficient 022,21 is at 
least a factor of 5 smaller than the corresponding one in tJi^, U4, Ud in any RDIs model. 
Hence, in improved models U22 and not Uim is the optimal quantity. 



5. 7. The critical exponent v 

We estimate the critical exponent v by using the data at p = 0.8, since in this case there 
are no leading corrections to scaling — the model is improved — the data have smaller 
errors, and we have results for larger lattices. We analyze the derivatives B!^ and U'i^ at 
fixed R^. 
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Figure 9. Estimates of U*^-^ obtained from fits of Uim to U*^ + bL The dashed 
hnes correspond to the estimate = 1.840(4) obtained from fits of U4, U22, Ud at 
p = 0.8. 



Since we have established that the RSIM at p = 0.8 is improved, the dominant 
scahng corrections are associated with the next-to-leading exponent uj2- Therefore, we 
fit R' to aL^/'' and to 

aL""" (1 + cL-"^) (52) 

or, more precisely. In i?' to In a + Mn L or to In a + Mn L + cL~'^^ . The exponent U2 is 
fixed to the FT value uj2 = 0.82(8). The results are reported in Fig. [10] as a function of 
Lmin. These fits are quite good, with x^/DOF < 1 already for relatively smal values of 
Linin. The dependence on lj2 is negligible when U2 varies within the range allowed by 
the FT estimate. The fits of R'^ to a simple power law are quite stable. The results for 
Lmin > 48 provide the estimate 

zy = 0.6835(10). (53) 

The exponent u can also be determined by analyzing the ratio x'/x? where x' is the 
derivative of x with respect to (3. This ratio also has the asymptotic behavior fl52|) . The 
corresponding results are in perfect agreement with those obtained by using R'^ and f/4. 

Collecting results we obtain the final estimate 

u = 0.683(2), (54) 

which takes into account the fits of R'^ and U^, with and without the scaling correction 
with exponent 002- In the determination of the estimate (JMl), we have implicitly assumed 
that the RSIM at p = 0.8 is exactly improved so that there are no leading scaling 
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Figure 10. Estimates of tlie critical exponent ly, obtained by fitting the MC data 
for the RSIM at p = 0.8 to a simple power law (above), i.e. to aL^^^ , and to Eq. ([5^ 
with UJ2 = 0.82. We consider i?^, U'^, R'^ and C/4 Some results are slightly shifted 
along the x-axis to make them visible. The dotted lines correspond to the final estimate 
V = 0.683(2). 



corrections. However, p* is only known approximately and thus some residual leading 
scaling corrections may still be present. To determine their relevance, we use the upper 
bound (H3l) and the MC results for R' aX p = 0.65. If the leading scaling corrections do 
not vanish, R' behaves as 

R! = aL'/" (1 + bR>,nL-^ + &KM2i^"'" + &fl',2ilv""^ + ■ ■ ■) • (55) 

In the following section we obtain the estimates 

= 0.65) = 0.60(15), bu^^nip = 0.65) = 0.40(15). (56) 

Bound (is]) gives then |6i?',ii(p = 0.80) | < 0.02 for both R'^ and U'^. We have thus 
repeated the analysis at p = 0.8 considering R'/{1 ± 0.02L^'^), with uj = 0.33. The 
results for the exponent u vary by ±0.0004, which is negligible with respect to the final 
error quoted in Eq. ( l54l) . 

The result ( 1541) is in agreement with and improves earlier MC estimates: u = 
0.683(3) (Ref. [T9]) and u = 0.6837(53) (Ref. [18]). It is also in agreement with the FT 
result tnj z/ = 0.678(10). 

To check universality, it is interesting to compute u directly, using the data at 
p = 0.65. We fitted \nR' to a + (l/u) InL + cL-"^ + dL'' taking u = 0.33(3) and 
e G [0.6, 0.9] (as before the last term takes into account corrections of order L^^'^ and 
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L~^^). We obtain v = 0.65(2), 0.67(2), 0.68(2) for L^i^ = 8, 12, 16. The somewhat large 
error is mainly due to the variation of the estimate with e. Thus, if scaling corrections 
are taken into account, universality is satisfied. 

5. 8. Improved estimators for the critical exponent v 

As we did in Sec. l5.6l for the phenomenological couplings, we wish now to define improved 
estimators of the critical exponent z/, i.e., quantities Bf-^^ such that = 0, see 

Eq. (!55|) . We will again combine data at p = 0.8 with data at p = 0.65. 

Let us consider a phenomenological coupling. In the following we choose IJd defined 
in Eq. f|T7j) because it has the least relative statistical errors among the combinations of 
f/4 and f/22. For generic values of p it has the asymptotic behavior 

Ud{p- V) = U: (1 + Q,n(p)^-" + ■■■). (57) 

Analyzing the data as described in Sees. 15.21 and [5.3[ we obtain the estimate Cd^uijp = 
0.65) = —0.16(2). Then, given a generic coupling R' with asymptotic behavior (|55l) . we 
have 

f^^^.d.M.,..-"....). (58, 

where Abn^^n = &i?Mi(p = 0.65) - 6r/,ii(p = 0.8). A fit of the MC data to Eq. ([58]) gives 
A^^n = 0.60(15) and Abu^^^u = 0.40(15). Bound implies \bR,^nip = 0.8) | < 
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\bR'.n{p = 0.65)1/30. Therefore, given the error bars on the previous results, we 
can identify A6r/ ii with = 0.65), obtaining the estimates aheady reported 

in Eq. ( 156|) . Finally we consider 

R^ = R'US ~ + (6^, 1, + aQ,n)i^-" + ■••]• (59) 

If we fix a = — the quantity is improved. Using the results obtained 
above, we find that 

Rlim^R'fiT With a5 = 4(l), (60) 
Kim = KUd'' with a, = 2.5(1.0) (61) 

are approximately improved quantities. As a check of the results, we perform the fits 
(1581) also for the improved observables. In both cases Abpi'^n is consistent with zero. 

At p = 0.8 these improved quantities give results which are substantially equivalent 
to those of the original quantities, as shown in Fig. [10], where we report the estimates 
corresponding to the central values of the exponents and a^. This is not unexpected, 
since the RSIM at p = 0.8 has suppressed leading scaling correction for any quantity. 

The use of the improved estimators is particularly convenient at p = 0.65. Indeed, 
as discussed in the previous Section, the direct analysis of R^ and U'^ gives estimates 
of u with a large error. In the improved quantities -R^ and im the leading scaling 
correction is absent and thus one should be able to determine more precise estimates of 
u and perform a more severe consistency check of universality. Since = 0, we fit 

the data to 

R^ = aL'/" (1 + cL-') (62) 

with e G [0.6,0.9]. The results of these fits are shown in Fig. [TT] (we report results 
corresponding to the two choices e = 0.66 ~ 2a; and e = 0.82 ?a u!2). Explicitly, we find 
u = 0.687(7), 0.684(7), 0.679(6), 0.680(8) from the analysis of R^^^ with L^^^ = 12, 16, 
24, 32. The error includes the statistical error, the error on a^, and the variation of the 
estimate with e. The results for U'^^^ are perfectly consistent. The estimates are three 
times more precise than those obtained by using the unimproved R^ and and are in 
good agreement with the more precise estimate u = 0.683(2) obtained from the data at 
p = 0.8. 



5. 9. Estimate of the critical exponent r] 

In order to estimate the critical exponent t] we analyze the magnetic susceptibility x- 
We have fitted the data at p = 0.8 to aL"^"^, to aLF'~^{l + cL~'^), taking e as a free 
parameter, and to aL'^~^{l + cL^'^^) with uj2 = 0.82,0.74,0.90 [we use, as before, the 
FT estimate uj2 = 0.82(8)]. Moreover, we simultaneously fit the data at p = 0.8 and 
p = 0.65 to aiL2-'?(l + ciL-°-S2) jf ^ ^ 9.8 and a2L2-''(l + caL-^-^S) H p = 0.65. The 
results are shown in Fig. [T2j They are fully consistent and provide the final estimate 

77 = 0.036(1). (63) 
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Figure 12. Estimates of tlie critical exponent rj, obtained by fitting the magnetic 
susceptibility x- All fits refer to the RSIM at p = 0.8, except the last one where we 
simultaneously fit data at p = 0.8 and at p = 0.65 (see text). Some results are slightly 
shifted along the x-axis to make them visible. The dotted lines correspond to the final 
estimate rj = 0.036(1). 



Again we should discuss the error due to possible residual leading scaling corrections at 
p = 0.8. A fit of the data at p = 0.65 to aL'^''^{l + cL~'^^) gives rj ^ 0.035. Even if we 
neglect the correction term proportional to , the result of the fit is compatible with 
the estimate (163|1 . This means that the amplitude of L""^ is quite small for p = 0.65 
and gives a correction of the order of the statistical error in Eq. (l63l) . The amplitude 
of at p = 0.8 is at least 30 times smaller, see Eq. (H3|) . Therefore, residual scaling 
corrections are negligible. 

The result fl63l) improves earlier estimates: rj = 0.035(2) and rj = 0.0374(45) from 
MC simulations of Refs. [I9] and [H] respectively. It is also close to, though not fully 
compatible with, the FT result t] = 0.030(3) [TT] . 

6. Finite-size scaling analysis of the random site-diluted Ising model at 
fixed P 

In this section we perform a different analysis using the results at /3 = /?run (the value 
of (3 at which we performed the MC simulation) for the RSIM at p = 0.8. Thus, the 
data we use here are different from those considered in the previous Section. We will 
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0.59403(31) [8] 
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Table 6. Results of the combined fits for the phenomenological couplings for several 
values of imin, the smallest lattice size used in the analyses. In the first set of fits, e is 
free, while in the second set we fixed e = ui2 = 0.82(8). The number in parentheses is 
the statistical error, while the number in brackets gives the variation of the estimate 
as UJ2 is varied within one error bar. 

combine them with thos^ obtained close to the critical point (0.28572 < (3 < 0.28578) 
in Ref. [19]. The statistics of the old results is comparable with that obtained here for 
p = 0.8; note however that we have now also results for the larger lattice size L = 192. 

6.1. Phenomenological couplings and critical temperature 

As discussed in Sec. I3.1[ close to the critical point a phenomenological coupling behaves 
as [see Eq. (ES])] 

R{(3, L) = R* + ai {(3 - (3^)1^^" + a2L-^ (64) 

In order to determine R* and jSc we have performed two different types of fit, always 
fixing u = 0.683(5) [results are essentially unchanged if we vary u in [0.678, 0.688], which 
is quite conservative given the estimate (l54l)]. In the first case, we simultaneously fit R^, 
t/4, and U22, keeping e as a free parameter. For the exponent e we obtain e = 0.95(20): as 
expected there is no indication of a correction-to-scaling term with exponent u ^ 0.33. 
In the second fit, we use the fact that the model is improved, set e = uj2, and use the 
FT estimate UJ2 = 0.82(8). The results of the fits for several values of L^in are reported 
in Table El They are quite stable with respect to Lmin and indeed results with L^i^ = 8 
are compatible with all those that correspond to larger values. If we take conservatively 
the final estimates from the results with Lmin = 24 and uj2 fixed, we have 

= 0.2857433(5), 
Rl = 0.5944(7), 
UI = 1.648(2), 

f/*2 = 0.1479(6). (65) 

% More precisely, the results of Ref. [19] consists in 5 data with L = 128 {Ng — 14000), 7 data with 
L = 64 and = 20000 and 3 data with L = 64 and iV, = 40000, 8 data with L ^ 32 and = 35000, 
and 4 data with L = 16 and Ng = 80000. Note that in Ref. 19J the derivatives R' were not determined, 
so that the analysis of v relies mainly on the present data. 
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Figure 13. Estimates of the critical exponent i^, obtained by simultaneous fits of 
and R'^, and of U4 and C/4 (see text). Some results are slightly shifted along the 
X-axis to make them visible. The dashed lines correspond to the estimate ly = 0.683(2) 
obtained in Sec. 15.71 



The estimate of Pc is compatible with that reported in Sec. I5.5[ Pc = 0.2857429(4). The 
estimate of i?| is essentially identical to that reported in Ref. [I9], = 0.5943(9), so 
that our analysis at fixed = 0.5943 corresponds indeed to fixing R^ = R'^. This is 
also confirmed by the results for and U22 that are compatible with the estimates of 
and U22 obtained in Sec. 15.21 (if we had performed analyses at fixed R^ 7^ R^ such an 
equality would not hold, see Sec. 13.21) . The estimates of and f/22 agree with previous 
MC estimates: U* = 1.650(9) and f/*2 = 0.1480(10) (Ref. ^LSj); = 1.653(20) and 
^72*2 = 0.145(7) (Ref. [H]). 



6.2. Estimates of u 

We now compute the critical exponent z/. It may be obtained by fitting R'{P = Pc, L) to 
aL^/^ This requires fixing Pc and this induces a somewhat large error. We have found 
more convenient to follow a different route, analyzing simultaneously R' (this gives v) 
and R (this essentially fixes P^. We performed two types of fits. First [fit (a)], we fit R 
and R' to 

R{p, L) =R* + aiip - pc)L^^'' + a2L-' 

+ asiP - PcfL^/" + a,iP - P,)L^'''-% 
R'iP, L) = aiL'/" + 2a^{p - pcjL^^" + aiL^'"-'. (66) 
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Figure 14. Estimates of the critical exponent 77, obtained in different fits (see text). 
Some results are slightly shifted along the x-axis to make them visible. The dashed 
lines correspond to the estimate 77 = 0.036(1) obtained in Sec. 15.91 



Here /3c and v are kept as free parameters, while e is fixed: e = LO2 = 0.82(8). In the 
second fit [fit (b)], we fit R and Bl to 

i?(/3, L) = i?* + ai(/3 - + a2L-^ 

ln[i?'(/3, L)] = 03 + ^ In L + a^ifi - fi^L^''' + a^L'^ (67) 

where again e = 102 = 0.82(8). The two fits give similar results, see Fig. [131 For instance, 
for Lmin = 24 and i? = i?^ we have v = 0.6814(17) and 0.6814(16) from fit (a) and (b), 
respectively (the reported errors are the sum of the statistical error and of the variation 
of the estimate as uj2 is varied within one error bar). If i? = f/4 we obtain analogously 
u = 0.6839(15) and 0.6840(16). Collecting results, this type of analysis gives the final 
estimate 

u = 0.6825(25), (68) 
which is fully compatible with Eq. fl5^ . 



6.3. Estimates of r) 

As in Ref. [19] we determine 77 from the critical behavior of Z = x/^^, which is more 
precise than x'- The relative error on x is 3.4 times larger than the relative error on Z. 
We perform two types of fits. First, we analyze Z as 

ln[Z(/3,L)] = a-r/ In L + 6Li/^(/3-/?J + cL-'^^ (69) 
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fixing u, Uii2, and (3c- The estimates are little sensitive to u and uj2 that are fixed to 
u = 0.683(5) and uj2 = 0.82(8). The dependence on is instead significant, of the 
order of the statistical error. We use Pc = 0.2857431(6) that combines the estimates 
determined in Sec. 15.51 and 16.11 Results are reported in Fig. [HI For Lmm = 24, 
Tj = 0.0364(9 + 2 + 8) where we quote the statistical error, the variation of the estimate 
with a; and p, and the change as jSc varies within one error bar. 

As done in Sec. 16.21 we can avoid using jSc by analyzing In Z together with a 
renormalized coupling R. We fit R to Eq. (Ell), again fixing z/ and e = UJ2- The results 
are reported in Fig. [HI For Lmm = 24 we obtain t] = 0.0367(15) and r] = 0.0367(13) 
using R^ and f/4, respectively. Collecting results this analysis provides the final estimate 

ri = 0.0365(15), (70) 

which agrees with the estimate r] = 0.036(1) obtained in the analyses at fixed R^. 

7. Finite-size scaling analysis of the randomly bond-diluted Ising model 

In this section we analyze the MC data of the RBIM at p = 0.55 and p = 0.7 at fixed 
R^ = 0.5943. They are reported in Tables H] and [5l respectively. We show that their 
critical behavior is fully consistent with that obtained for the RSIM in Sees. [5] and [6] 

7.1. Phenomenological couplings 

We first discuss the FSS behavior of the quartic cumulants. In Fig. [6] we have already 
shown the MC results U22{L) for the RBIM at p = 0.55 and p = 0.7 versus L~'^ with 
uj = 0.33. Their large-L behavior is perfectly consistent with that observed in the RSIM, 
all data converging to U22 = 0.148(1) as L increases. A more quantitative check can 
be performed by fitting the MC results for Uim- Indeed, as discussed in Sec. 15. 6[ the 
leading scaling corrections in f/im are at least a factor of 10 and a factor of 20 smaller 
that those in U22 and Ud, respectively. Thus, for any generic p, we should be able to 
obtain estimates of If*^ that are more precise than those of U22 and U^. In particular, 
we expect errors comparable to that of the RSIM result U*^ = 1.840(4). We fit f/im 
to U*^ + cL~^ with e = 0.66 and 0.82 (we remind the reader that the leading scaling 
corrections to are proportional to and L~^'^). Results are reported in Fig. [9] 
They are fully compatible with those obtained in the RSIM at p = 0.8, confirming that 
the RBIM belongs to the same universality class of the RSIM. 

A direct analysis of U22 and Ud for the RBIM at p = 0.7 gives results with 
large errors (see the corresponding analysis for the RSIM aX p = 0.65 reported in 
Sec. 15. 3p . Universality is verified, though with limited precision. More precise results 
are obtained for the RBIM at p = 0.55, since for this value of p the RBIM turns out to 
be approximately improved. To determine p* for the RBIM we follow the same strategy 
employed in Sec. 15.41 We determine the correction-to-scahng amplitude 022,11 obtaining 
C22,ii(p = 0.55) = —0.01(2) and 022,11 (p = 0.7) = —0.17(2). Then, we assume that for 
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Figure 15. Estimates of U22, U*^-^, and [7j obtained by fitting U22, Uim, and Ud for 
the RBIM at p = 0.55, to U* + cL''^^ with W2 = 0.82. The dashed hnes correspond to 
the estimates of tj* obtained from the analyses of the data of the RSIM a,t p — 0.8. 
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1.07(15). 



(71) 



(72) 



(73) 



for the RBIM. Again, in setting the error we have not considered the error on the hnear 
interpolation (172|) . 

Since the RBIM at j9 = 0.55 is approximately improved, we can fit [722, t^im, and 
Ud to U* + cL~^'^ . In Fig. [15] we show the results for uj2 = 0.82. They are very stable 
and in perfect agreement with the results obtained from the FSS analysis of the RSIM 
at p = 0.8. Indeed, we obtain 



where the error takes into account the uncertainty on uj2. They must be compared 
with the estimates obtained from the FSS analysis of the RSIM: U*^ = 1.840(4), 
U22 = 0.148(1), and = 1.500(1). These results provide strong evidence for 
universality between the RSIM and the RBIM. 

7.2. Critical temperatures 

Let us now estimate f3c from the estimates of l3f. For p = 0.55, since the model is 
approximately improved, we can fit the data of /?/ to pc + cL-^^"-^''. We obtain = 
0.4322895(15), which is compatible with the MC results of Ref. UM, Pc = 0.43225(10). 
The analysis of the 19th-order high-temperature expansion of x reported in Ref. [8j gave 
the estimate p^ = 0.43253(12). 

For p = 0.7 we must take into account the leading scaling corrections, i.e. fit Pf 
to Pc + cL~^, where e G [cu + l/h',U2 + 1/z^]. We obtain Pc = 0.326707(2), which agrees 
with the MC estimate pc = 0.32670(5) of Ref. [2T]. 

7.3. Critical exponents 

Let us now consider the critical exponents. Since the RBIM at p = 0.55 is approximately 
improved, we perform the same analysis as done for the RSIM at p = 0.8, see Sec. 15.71 
In Fig. [16] we show the results of several fits of the data of the derivative of the 
phenomenological couplings at p = 0.55. The estimates obtained by using R'^ and R'^ ^^ 
are substantially equivalent, as expected because the RBIM at p = 0.55 is approximately 
improved. The estimates of u shown in Fig. [16] are fully consistent with the estimate 
v = 0.683(2) obtained from the FSS analysis of the RSIM at p = 0.8. Similar conclusions 
hold for the critical exponent rj. In Fig. [17] we show the results of several fits analogous 



U*^ = 1.842(4), U;^ = 0.148(1), U*, = 1.501(2), 



(74) 
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Figure 16. Estimates of the critical exponent ly, obtained by fitting the data of the 
RBIM at p = 0.55 (some resuhs are slightly shifted along the a:-axis to make them 
visible). The dashed lines correspond to the estimate v = 0.683(2) obtained in the 
RSIM at p = 0.8. 

to those discussed in Sec. 15. 9[ Again universality is well satisfied: the estimates of r] are 
compatible with the RSIM result rj = 0.036(1). 

In the case of the RBIM at p = 0.7, analyses of unimproved quantities give estimates 
of u with large errors. We therefore only consider the improved quantities. In Fig.fTSlwe 
show the results of fits of R'^^^ for the RBIM at p = 0.7. They are again substantially 
consistent with the estimate u = 0.683(2) obtained from the FSS analysis of the RSIM 
at p = 0.8. The MC estimates of U'^^i-^ are less precise, but again consistent with 
universality. 
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Appendix A. Field-theory estimate of u)2 

In this appendix we estimate uj2 by a reanalysis of its FT six-loop expansion in the 
massive zero-momentum scheme. In Ref. [11] we performed a direct analysis of the 
stabihty matrix at the RDIs fixed point (FP), obtaining UJ2 = 0.8(2). Here we shall use 
the method discussed in Ref. [19] — it consists in an expansion around the unstable Ising 
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Figure 17. Estimates of tlie critical exponent 77, obtained by fitting the data of the 
RBIM at p = 0.55. The dashes line correspond to the estimate 77 — 0.036(f) obtained 
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FP — which allows us to estimate accurately the difference 002 — oois, where ujis is the 
leading correction-to-scaling exponent in the standard Ising model. Since uis is known 
quite precisely, this method allows us to obtain a precise result for lj2. 

In the FT approach one starts from Hamiltonian (jlj), determining perturbative 
expansions in powers of the renormalized couplings u and v. We normalize them so that 
u ^ Uq and f ~ f at tree level (these are the normalizations used in Ref. [19]; they 
differ from those of Ref. [H]). As discussed in Ref. [9], it is more convenient to introduce 
new variables y = u + v and z = —u. The Ising FP is located at = and y] = g^^, 
where p] g^ = 23.56(2). The RDIs FP is located at y* = 24.7(2) and z* = 18.6(3) (we 
use the MC results of Ref. pLSj since the FT estimates y* = 24.8(6) and z* = 14(2) are 
less precise). 

The expansion around the Ising FP can be performed along the Ising-to-RDIs RG 
trajectory [9], which is obtained as the limit Zq 0^ {uq 0^) of the RG trajectories 
in the z, y plane, see Fig. [H An effective parametrization of the curve is given by the 
first few terms of its expansion around z = 0, which is given by 

y-yi = T{z) = C2Z^ + c^z^ + ■ ■ ■ (A.l) 

where [9i| C2 = 0.0033(1) and C3 = 1(2) x 10~^. The fact that y — yi is of order z^ is the 
main reason why the variable y was introduced and is due to the identity [9] 



du 



11=0 



013., 



w=0 



= 0, (A.2) 

M=0 
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Figure 18. Results for the critical exponent ly, obtained by fitting the data of the 
RBIM at p = 0.7 (some results are slightly shifted along the x-axis to make them 
visible). The dashed lines correspond to the estimate ly = 0.683(2) obtained in the 
RSIM at p = 0.8. 



which corresponds to 



dz 



= 0. (A.3) 

z=0 

Performing the variable change y = g + T{z) in the double expansion of a generic 
quantity f{y, z) in powers of y and z, we obtain 

f{g, z) = f{g + T{z), z) = ^ e,{g)z\ (A.4) 

i 

The coefficients must be evaluated at g = g^^. This is done by resumming their 
perturbative expansions as discussed in Ref. [30j: in particular, we exploit Borel 
summability and the knowledge of the large-order behavior at the Ising FP. 

In Ref. [19] this approach was applied to the standard critical exponents. Here we 
extend these calculations to the next-to- leading scaling-correction exponent 002- For this 
purpose, we consider the stability matrix 

Q = {df3y/dy, df3y/dz; df3Jdy, d(3Jdz). (A.5) 



Each entry has an expansion of the form ( ]A.4[) . with coefficients ei{g) that are resummed 
as discussed before. Then, the matrix Q is diagonalized, obtaining the expansion of its 
eigenvalues in powers of z up to 0{z^). The corresponding coefficients for the smallest 
eigenvalue are quite large, so that this method does not provide an estimate of u which 
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is more precise than that obtained in Ref. [TI], i.e. u = 0.25(10). On the other hand, 
the expansion coefficients for uj2 are quite small. To order we obtain 

UJ2 - (^Is = '^CiZ\ (A. 6) 

i 

ci = 0, C2 = 5(2) X 10~^ C3 = -2(6) x 10-^ (A.7) 

where Ci = exactly [this is a consequence of relation flA.Sp ]. and the errors on the 
coefficients C2 and C3 are due to the resummation of the corresponding series evaluated 
at Qj^. Expansion (1A.6P is evaluated at z = 18.6(2), obtaining 



cu2-cuis = 0.00(5), (A.8) 

where the error takes into account all possible sources of uncertainties: the error on 
the coefficients, the truncation of the series in powers of z, and the uncertainty on the 
coordinates of the RDIs FP. Then, using the estimate uis = 0.82(3), which takes into 
account the results of Refs. [161 EH [39] obtained by various approaches, we arrive at the 
estimate 

UJ2 = 0.82(8), (A.9) 
which improves the result UJ2 = 0.8(2) obtained in Refs. pT| [9]. 

Appendix B. Bias corrections 

In this section we discuss the problem of the bias correction needed in the calculations 
of disorder averages of combinations of thermal averages. As already emphasized in 
Ref. [38], this is a crucial step in high-precision MC studies of random systems. 

To discuss it in full generahty, let us indicate with S the state space corresponding 
to the variables a and with R that corresponding to the dilution variables. Then, we 
consider a probability function 7r((j; p) on S depending parametrically on p (vr = e^^^/Z 
in the specific calculation) and a probability p{p) on R. Averages over TT{cr] p) are 
indicated as (■), or with (■)p when we wish to specify the value p used in the calculation. 
Moreover, we assume R to have a finite number K of elements {K = 2^ and K = 2'^^ 
in the RSIM and in the RBIM, respectively). We wish to compute averages of functions 
A{a,p). We ffist discuss the calculation of 



a=(A)" = J]p(p) 



^n{a;p)A{(x, p) 



(B.l) 



A numerical strategy to compute On could be the following. Extract A^^ independent 
disorder configurations pa, a = 1, . . . , Ng, with probability p{p) and then, for each p^, 
extract Nm independent configurations aa^a, cl = 1, • • • , A^m, with probability n^a; pa)- 
Then, define the sample average 

^ Nm 



a=l 
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A possible estimator of C„ could be 

1 

* a=l 

The question is whether O^^ converges to On defined in Eq. ( IB.ip as A^^ — oo at fixed 

Nm. 

To answer this question, let N{p) be the number of pa such that pa = p- Eq. (]B.3p 
can thus be rewritten as 

pGiJ a=l 

where the second sum extends over the N{p) terms that appear in Eq. (1B.3P such that 
Pa = p, i.e., which correspond to the same disorder configuration. As Ng — > oo, N{p) 
converges to Nsp{p) with probability 1; thus, as A^^ ^ oo 

Nip) Nip) 

Thus, for A^s ^ oo we have 

oT = Y.p^p)m;)^ = WFi. (B.6) 

p 

Eq. (1B.6P relies only on the limit Ng oo and is valid for any A^^- If n = 1 we obtain 

\a=l / p 

SO that Of* converges to Oi irrespective of A^: one could even take = 10 
This result could have been guessed directly from Eq. (IB.ip . since 

Oi = J2p{p)n{a-p)A{a,p). (B.8) 
p 

A correct sampling is obtained by determining each time a new a and p with combined 
probability p{p)'K{a\ p). Let us now consider n = 2. We have 

^ / N„, Nm ' 



ATS 



™ \a=l 6=1 



-i^ [Ar^(A^^ - 1) {A)l + Nm {A')^ 

m 



In practice, one could fix Nm in such a way to minimize the variance 



varOf = {^/Ns)[{A)^ - {{A)f + {{A^) - {A)^)/Nm]. 
For X st fixed /3 this minimization can be done explicitly. Indeed, the variance can be related to 
174 and C/22: (err[x])^/x^ = [U22 + {Ua - \)/Nm\/Ns. Thus, at the critical point (err[x])^/x^ = 
(0.148 + 0.648/A'm) /Ns- If the work for each disorder configuration is proportional to A^thorm + Nm and 
-^thcrm — 300, it is casy to verify that the optimal Nm (the one that gives the smallest errors at fixed 
computational work) corresponds to Nm = 53. 
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Thus, we obtain for A^^ — > oo at fixed N^n 

Of ^02 + -^ (W) - 'W) ■ (B. 10) 



The second term is what is caUed the bias. Since in the simulations Nm is finite and not 
too large, this term may give rise to systematic deviations larger than statistical errors. 
It is therefore important to correct the estimator in such a way to eliminate the bias. 

For this purpose we divide the configurations in n bunches and define the 
sample average over bunch i of length N^/ n: 

iNm/n 

[A]l/n,i,p„ = ^ J2 M(^a,a,Pa). (B.ll) 

a=l+{i-l)Nm/n 

A new estimator of On is 

1 

^unbiased ^ 5]][^]l/n,l,p„ [^]l/n,2,p„ " " " [^]l/n,n,p„- (B-12) 



s 



a=l 

Let us verify that this estimator is unbiased. By repeating the arguments presented 
above, for Ng ^ 00 the estimator (^^^ibiased converges to 



^unbia^ed ^ (^[A]yn,M]l/n,2 ■ ■ ■ [A]l/n,n>. (B.13) 

Because the configurations are assumed to be independent, we have 

([A]i/„,i[A]i/„,2 ■ ■ ■ = ([A]i/„,2> ■ ■ ■ {[A]i/n,n) = {Af . (B.14) 

Thus, for any n, irrespective of A^^ (one could even take = n), C)^ribiased converges 
to On as Ng — > 00. 

The considerations reported above can be trivially extended to disorder averages 
of products of sample averages. Thus, in order to compute {A){B), we use 

1 

^ ([^]l/2,l,p„[-B]l/2,2,p„ + [B]i/2,1,pAA]i/2,2,p^) , (B.15) 



while for {A){B){C) we use 



SW ^ ^ ^Z^'^'"" ^Z^'^'"" ^Z^'^'"" ^ ^ permutations} . (B. 16) 

* a=l 

In the case of n terms, we divide the estimates into n parts and then consider all 
the n\ permutations. 

In this paper we extensively use the reweighting technique that requires the 
computation of averages of the form 



Ra,b={^], (B.17) 



where the disorder average should be done after computing the ratio. Indeed, given a 
MC run at inverse temperature /5 the mean value of an observable O at f3 -\- A/5 is given 

by 

(0),+A, = {Oe-^n,/{e-^np- (B.18) 
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A.B 



^ a=l \ 



1/2,1,Pq 



1 



l/2,2,pc 



40 



(B.19) 



l/2,l,Pc, 



(B.20) 



where {■} should be defined so that 

-I- — 

Bl (B) 

for Ng ^ oo and any (the suffix 1/2, i has the same meaning as before). We have 
not been able to define an estimator with this property. We thus use a biased estimator, 
with a bias of order N~'^. Consider 



1 



l^fB{aa,p)-{B), 
N„- ^ 



(B.21) 



Assuming large, we can expand the term in brackets, keeping the first nonvanishing 
term: 



1 



1 



1 



ab 



N, 



1 {B')^{Br 

{B)l 



Thus, if we define 
1 

B 

which is such that 



[B] 



1 m-iBi^ 



1 



{B) 



;i + o(Ar„2)), 



+ ■ ■ 



(B.22) 



(B.23) 



(B.24) 



the estimator ( IB.IQP converges to Rab with corrections of order . Thus, when using 
the reweighting technique, Nm is crucial and cannot be too small. 

In this paper we also compute derivatives of different observables with respect to 
j3. They can be related to connected correlation function as 



d{0) 
dp 



-{on)-{o){n). 



(B.25) 



Therefore, if we apply the reweighting technique, we need to compute terms of the form 

(B.26) 



Rab. 



c 



{A){B) 



with A = Oe~^'^'^, B = TYe"'^^^, and C = e'^'^'^. A possible estimator (this is the one 
that is used in the paper) is 

Ns 



AJ Z-^ 



a=l 



[^]l/4,l,p„ [5] 1/4,2, 



Pa 



1 



^ J l/4,3,p<, I J 1 



1 



permutations 



L/4,4,pc 



.(B.27) 
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The formulae we have derived above rely on two assumptions: (i) different configurations 
obtained with the same disorder pa are uncorrelated; (ii) configurations a are extracted 
with probability 7r{a,p). None of these two hypothesis is exactly verified in practical 
calculations. Hypothesis (i) is violated because MC simulations usually provide 
correlated sequences of configurations. Correlations change some of the conclusions 
presented above. First, the estimator flB.12|) is no longer unbiased, except in the case 
n = 1. To explain this point, let us consider the specific case n = 2. In the presence of 
correlations 

(AK,,p,)A(a,„,p,))^^ = {A)l + {Y^T^^A)C{\a-b\), (B.28) 
where vslt A = (A)'^ — (A)'^ and C{t) is the autocorrelation function. Then 



/.\2 , 4varp^v4 



^ aC{a) + ^ {Nm - a)C(a) 

a=l a=iV™ /2+1 



(B.29) 



If C{t) decays fast enough the term in brackets is finite for Nm oo. If we further 
assume C{t) = exp(— t/r), we find that the bias is of order (r/A^^)^. Thus, for n > 2 
the estimator flB.12p is biased. Therefore, it is no longer possible to take Nm at will. 
To avoid the bias, A^^ should be significantly larger than the autocorrelation time. 
Analogously, Eq. (IB. 231) is correct only in the absence of correlations. Otherwise one 
should define 

[b}-W]V-^ — m^)- 

where tb, mt is the integrated autocorrelation time associated with B. In the present 
work the MC algorithm is very efficient, so that our measurements should be nearly 
independent. Thus, we have always used Eq. (lB.23p . 

The second assumption that we have made is that configurations aa^a are extracted 
with probability 7i{a,pa). In MC calculations this is never the case: configurations 
are obtained by using a Markov chain that starts from a nonequilibrium configuration. 
Therefore aa,a are generated by using a distribution that converges to 7r(cr, pa) only 
asyniptotically. This is the so-called inizialization bias. This bias, which is of order 
cannot be avoided. However, by discarding a sufficiently large number of initial 
configurations one can reduce it arbitrarily. Note that, as A^^ is increased, either or 
the number of discarded configurations should be increased too. 



* More precisely, if we define 



then Nm{{[A]Nti^) MC ~ (^}) converges to a nonvanishing constant K as ~^ oo (here {■)mc is an 
average over all Markov chains that start at i = from an arbitrary configuration). The constant K 
decreases as TVth increases, with a rate controlled by the exponential autocorrelation time. 
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Figure Bl. Ui and R'^ for a run of the RSIM model at L = 64, p = 0.8, TV, = 60000, 
/3run — 0.285742, reweighted at = 0.5943 as a function of l/N,n- In the insets we 
show the results for the bias-corrected estimates versus l/N"^. We report data with 
iV™ = 40,80,120, ...,400. 



Of course, the practicaUy interesting question is whether, with the values of Nm 
and Ng used in the MC simulations, bias corrections are relevant or not. For this 
purpose, in Fig. IBll we compare, for a specific run of the RSIM at p = 0.8, the value 
of the derivatives R'^, with and without bias correction, as a function of Nm] data 
obtained for /? = 0.2857420 are reweighted to have = 0.5943 [this corresponds to 
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(3f = 0.2857478(18)]. The average without bias correction has been estimated by using 



(B.31) 



while the bias-corrected estimate corresponds to Eq. (]B.27p . From the figure, we see 
that the biased estimate shows a systematic drift, which, as expected, is hnear in 1/Nm- 
The bias-corrected estimate is instead essentially fiat; deviations are observed only for 
f/4 for N.m ^ 100. As expected, they apparently decrease as Note, that even 

for Nm = 400 the biased estimate differs within error bars from the bias-corrected one. 
Thus, with the precision of our calculations, the bias correction is essential for both R'^ 
and t/4. 

We would like to point out that the expressions obtained here are somewhat different 
from those proposed in Ref. [3H]. Only for {A){B) are they identical. We report here 
the expressions for Rab and Rabc- 
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■A,B 
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1/2,1 



[A] 
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(B.32) 



4[C]^/2,l '*L'-'Jl/2,2 

The idea behind these formulas is the following. Consider O which is an arbitrary 

function of thermal averages. To compute O, consider an arbitrary estimator O"^^^ of O. 

For — >■ cx) at fixed Nm, O'^^^ converges to 
a 



+ 



+ Oil/Ni) 



A better estimator, without the l/N^ correction, is obtained by considering 

1 



est,unb 



20-* - ^or/2,1 - loi%,2- 



(B.33) 



(B.34) 



Here O^^ is determined by using all Nm measures, while C^^2 1 and Of^^^ 2 computed 
by using the first half and the second half of the measures, respectively. It is easy to 
see, by substituting the behavior (IB. 331) in Eq. ( IB. 341) . that the new estimator has no 
corrections of order 1/Nm- 
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